Models Update
Jina Park, Phuc Nguyen
12/5/2017
Progress Made
Made maps and graphs of our models
Altered our models to get them to converge.
Group Roles
- We worked together to make the maps/graphs and fixed the models.
Research Question 1
Predicting hotspots of specific crimes and how they changed over time and space.
Are certain types of crimes more prevalent in certain areas than others?
How have these hotspots changed over time? If they have changed over time, how has it changed relative to changes in income, demographics, age, etc?
Research Question 2
Examining police traffic stops and subsequent arrests to determine whether certain ethnic minorities are stopped/arrested more than others.
Are certain ethnic groups more likely to be arrested/stopped in certain areas than others? For instance, are blacks/hispanics more likely to be stopped/arrested in predominantly white neighborhoods?
We would like to examine the relationship between traffic stops/arrests and the actual crime rates of different ethnic groups. For instance, are blacks/hispanics getting stopped/arrested at a disproportionately greater rate than they are committing crimes?
Load Libraries
library(tidyr)
library(rjags)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(MacBayes)
library(choroplethrZip)
library(rgdal)Map cleanup
cleanup <- theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_rect(fill = 'white', colour = 'white'),
axis.line = element_line(colour = "white"),
axis.ticks=element_blank(), axis.text.x=element_blank(),
axis.text.y=element_blank())Clean data
#Load data
Crime = read.csv('Crime_sample_2010_to_2017.csv')
Crime = Crime %>%
mutate(Date = as.POSIXct(strptime(Date, "%Y-%m-%d %H:%M:%S"))) %>%
filter(!is.na(Latitude))## Warning in strptime(Date, "%Y-%m-%d %H:%M:%S"): unknown timezone 'zone/tz/
## 2017c.1.0/zoneinfo/America/Chicago'Project observations to ZIP codes boundaries
#load shapefile
zip_shp <- readOGR(dsn="Boundaries - ZIP Codes", layer="geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a")## OGR data source with driver: ESRI Shapefile
## Source: "Boundaries - ZIP Codes", layer: "geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a"
## with 61 features
## It has 4 fieldstemp <- zip_shp
zip_df <- fortify(temp, region = "zip")
#needs to reassign CRS for shapefile
new_CRS <- CRS("+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
zip_shp <- sp::spTransform(zip_shp, new_CRS)
#summarize crime count per year
years <- unique(Crime$Year)
Hotspot_Crime <- data.frame()
for(i in 1:length(years)){
year <- years[i]
Current_Year <- Crime %>% filter(Year == year)
locations <- with(Current_Year, as.data.frame(cbind(Longitude, Latitude)))
coordinates(locations) <- ~Longitude+Latitude
proj4string(locations) <- CRS("+init=epsg:4326")
#prepping the data: project points to polygons
proj4string(zip_shp) <- new_CRS
proj4string(locations) <- new_CRS
by_zip <- over(locations, zip_shp)
by_zip <- by_zip %>%
group_by(zip) %>%
dplyr::summarise(total=n()) %>%
filter(!is.na(zip)) %>%
mutate(id = as.character(zip))
total_map <- left_join(zip_df, by_zip)
total_map <- total_map %>% mutate(Year = year)
Hotspot_Crime <- rbind(Hotspot_Crime, total_map)
}Add demographic information to Model Data
data("df_zip_demographics")
Demographics <- df_zip_demographics %>% mutate(zip = region)
Hotspot_Crime <- inner_join(Hotspot_Crime, Demographics, by="zip") ## Warning: Column `zip` joining factor and character vector, coercing into
## character vectorModel_Data <- Hotspot_Crime %>%
select(-c(lat, long, order, piece, group, hole)) %>% distinctGet populations by zipcode and year
Crime_count <- Model_Data %>%
select(Year, zip, total) %>%
spread(key=Year, value=total)
rownames(Crime_count) <- Crime_count$zip
Crime_count <- Crime_count %>% select(-zip)
Crime_count[is.na(Crime_count)] <- 0
head(Crime_count)## 2010 2011 2012 2013 2014 2015 2016 2017
## 60601 8 6 5 8 7 6 7 11
## 60602 4 8 1 4 5 10 8 9
## 60603 1 5 3 2 1 2 9 4
## 60604 1 6 3 1 4 3 1 3
## 60605 10 7 13 6 9 10 10 8
## 60606 2 6 4 3 1 4 3 2Population <- Model_Data %>%
select(total_population, zip, Year) %>%
spread(key=Year, value=total_population) %>%
select('2010', zip)
rownames(Population) <- Population$zip
Population <- Population %>% select(-zip)
Population[is.na(Population)] <- 0
head(Population)## 2010
## 60601 9732
## 60602 1463
## 60603 882
## 60604 415
## 60605 24972
## 60606 2789Population.vector <- Population$`2010`
time <- sort(unique(Model_Data$Year))
percent_nonwhite <- 100 - Model_Data$percent_whiteSimulation
Research Question 1
Plot of \(log(y_{ij}/n_{i})\) vs year
Crime_count_log_calc <- log(Crime_count/Population.vector)
Crime_count_log <- Crime_count_log_calc %>%
rownames_to_column() %>%
gather(key=year, value=crime_count_log, -rowname) %>%
mutate(zip=rowname) %>%
select(zip, year, crime_count_log) %>%
mutate(year=as.numeric(year))
head(Crime_count_log)## zip year crime_count_log
## 1 60601 2010 -7.103733
## 2 60602 2010 -5.901950
## 3 60603 2010 -6.782192
## 4 60604 2010 -6.028279
## 5 60605 2010 -7.822925
## 6 60606 2010 -7.240291ggplot(Crime_count_log, aes(x=year, y=crime_count_log, colour=zip)) + geom_point() + geom_line() 
Plot of \(log(y_{ij}/n_{i})\) vs per capita income
Model_Data_1 <- Model_Data %>%
filter(Year == 2010) %>%
mutate(percent_nonwhite = 100 - percent_white) %>%
select(zip, per_capita_income, percent_nonwhite) %>%
right_join(Crime_count_log, by="zip")
head(Model_Data_1)## zip per_capita_income percent_nonwhite year crime_count_log
## 1 60601 89136 37 2010 -7.103733
## 2 60602 66748 38 2010 -5.901950
## 3 60603 73239 39 2010 -6.782192
## 4 60604 135807 20 2010 -6.028279
## 5 60605 58338 42 2010 -7.822925
## 6 60606 83005 25 2010 -7.240291ggplot(Model_Data_1, aes(x=log(per_capita_income), y=crime_count_log, colour=percent_nonwhite)) + geom_point() + labs(x="Log of Per Capita Income", y="Log of Crime Rate", title="Relationship between Per Capita Income and Log of Crime Rate")
Plot of \(log(y_{ij}/n_{i})\) vs nonwhite percentage
ggplot(Model_Data_1, aes(x=log(percent_nonwhite), y=crime_count_log)) + geom_point() + labs(x="Log of Nonwhite Percentage", y="Log of Crime Rate", title="Relationship between Nonwhite Percentage and Log of Crime Rate")
Model 0
Data
\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \beta * time_{j}\]
Priors
\[\beta \sim N(0, 1000^2)\]
i: zipcodes in Chicago (58 total in our data)
j: years (from 2010 to 2017)
y[i, j]: observed count of violent crimes in Chicago by year and zipcode
p[i, j]: probability of being a victim in location i and time j
n[i]: the population in zipcode i
\(\beta\): time trend
time[j]: the time in years
#Data
Crime_count <- Model_Data %>%
group_by(Year) %>%
dplyr::summarise(total.stop = sum(total))
head(Crime_count)## # A tibble: 6 x 2
## Year total.stop
## <int> <int>
## 1 2010 1474
## 2 2011 1448
## 3 2012 1340
## 4 2013 1258
## 5 2014 1090
## 6 2015 1117#Specify the model
crime_model_0 <- "model{
#Data
for(j in 1:8){
y[j] ~ dbin(p[j], n)
logit(p[j]) = beta*time[j]
}
#Priors
beta ~ dnorm(0, 1/100^2)
}"
#set up an algorithm to simulate the posterior by
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_0 <- jags.model(textConnection(crime_model_0), data=list(y=Crime_count$total.stop, n=sum(Population.vector), time=Crime_count$Year),inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 8
## Unobserved stochastic nodes: 1
## Total graph size: 40
##
## Initializing model#simulate a sample from the posterior
#note that we specify both mu and tau variables
crime_sim_0 <- coda.samples(crime_jags_0, variable.names=c("p", "beta"), n.iter=50000)
#store the samples in a data frame:
crime_data_0 <- data.frame(crime_sim_0[[1]])Running Mean Plot
running_mean_plot(crime_data_0$beta)
Parameter Inferences
ggplot(crime_data_0, aes(x=beta)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Effect of Time", x="Time Effect", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.1.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2010", x="Probability of Crime", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.2.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2011", x="Probability of Crime", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.3.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2012", x="Probability of Crime", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.4.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2013", x="Probability of Crime", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.5.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2014", x="Probability of Crime", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.6.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2015", x="Probability of Crime", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.7.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2016", x="Probability of Crime", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.8.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2017", x="Probability of Crime", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Confidence Intervals and Posterior Means of Crime Model 0 Parameters
crime_model_0_CI <- matrix(nrow=9, ncol=2)
for (i in 1:9) {
crime_model_0_CI[i,] <- quantile(crime_data_0[[i]], c(0.025, 0.975))
}
crime_model_0_mean <- matrix(nrow=9, ncol=1)
for (i in 1:9) {
crime_model_0_mean[i, ] <- mean(crime_data_0[[i]])
}
crime_model_0_table <- cbind(crime_model_0_CI, crime_model_0_mean)
rownames(crime_model_0_table) <- names(crime_data_0)
colnames(crime_model_0_table) <- c("2.5%", "97.5%", "Mean")
as.table(crime_model_0_table)## 2.5% 97.5% Mean
## beta -0.0038479706 -0.0038283390 -0.0038380756
## p.1. 0.0004373144 0.0004549075 0.0004461179
## p.2. 0.0004356356 0.0004531701 0.0004444097
## p.3. 0.0004339632 0.0004514393 0.0004427081
## p.4. 0.0004322972 0.0004497151 0.0004410130
## p.5. 0.0004306377 0.0004479975 0.0004393244
## p.6. 0.0004289845 0.0004462865 0.0004376422
## p.7. 0.0004273377 0.0004445820 0.0004359665
## p.8. 0.0004256971 0.0004428839 0.0004342971Model 1
Notation
i: zipcodes in Chicago (58 total in our data)
j: years (from 2010 to 2017)
y[i, j]: observed count of violent crimes in Chicago by year and zipcode
p[i, j]: probability of being a victim in location i and time j
n[i, j]: the population in zipcode i in time j.
time[j]: the time in years
\(\delta\)[i]: the spatial temporal interaction term
Specify the Model
We are modeling the count of crime incidents by year and zipcodes using a binomial distribution, in which the parameters are the probability of being a victim and the total population in location i and year j.
We modeled the probability of being a victim using a logit function of the grand crime rate of Chicago and the inclusion of several random effects. It is also a function of time and its interaction with location.
Data
\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \delta_{i} * time_{j}\]
Priors
\[\delta_{i} \sim N(0, prec.delta)\]
Hyperpriors
\[prec.delta \sim Gamma(0.5, 0.0005)\]
#Data
Crime_count <- Model_Data %>%
select(Year, zip, total) %>%
spread(key=Year, value=total)
rownames(Crime_count) <- Crime_count$zip
Crime_count <- Crime_count %>% select(-zip)
Crime_count[is.na(Crime_count)] <- 0
head(Crime_count)## 2010 2011 2012 2013 2014 2015 2016 2017
## 60601 8 6 5 8 7 6 7 11
## 60602 4 8 1 4 5 10 8 9
## 60603 1 5 3 2 1 2 9 4
## 60604 1 6 3 1 4 3 1 3
## 60605 10 7 13 6 9 10 10 8
## 60606 2 6 4 3 1 4 3 2#Specify the model
crime_model_1 <- "model{
#Data
for(i in 1:58) {
for(j in 1:8){
y[i,j] ~ dbin(p[i,j], n[i])
logit(p[i,j]) = delta[i]*time[j]
}
}
#Priors
for(i in 1:58){
delta[i] ~ dnorm(0, prec.delta)
}
#Hyperpriors
prec.delta ~ dgamma(0.5, 0.0005)
}"
#set up an algorithm to simulate the posterior by
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_1 <- jags.model(textConnection(crime_model_1), data=list(y=Crime_count, n=Population.vector, time=time),
inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 464
## Unobserved stochastic nodes: 59
## Total graph size: 1520
##
## Initializing model #simulate a sample from the posterior
#note that we specify both mu and tau variables
crime_sim_1 <- coda.samples(crime_jags_1, variable.names=c("p", "delta"), n.iter=60000)
#store the samples in a data frame:
crime_data_1 <- data.frame(crime_sim_1[[1]])Parameter Interpretation
ggplot(crime_data_1, aes(x=delta.1.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Crime Differential Growth in 2010", x="Crime Differential Growth", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_1, aes(x=p.1.1.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2010 in Zipcode 1", x="Probability of Crime", y="Density")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Confidence Intervals and Posterior Means of Crime Model 1 Parameters
crime_model_1_CI <- matrix(nrow=length(names(crime_data_1)), ncol=2)
for (i in 1:length(names(crime_data_1))) {
crime_model_1_CI[i, ] <- quantile(crime_data_1[[i]], c(0.025, 0.975))
}
crime_model_1_mean <- matrix(nrow=length(names(crime_data_1)), ncol=1)
for (i in 1:length(names(crime_data_1))) {
crime_model_1_mean[i, ] <- mean(crime_data_1[[i]])
}
crime_model_1_table <- cbind(crime_model_1_CI, crime_model_1_mean)
rownames(crime_model_1_table) <- names(crime_data_1)
colnames(crime_model_1_table) <- c("2.5%", "97.5%", "Mean")
as.table(crime_model_1_table)## 2.5% 97.5% Mean
## delta.1. -3.711729e-03 -3.456400e-03 -3.580096e-03
## delta.2. -2.867138e-03 -2.586290e-03 -2.722396e-03
## delta.3. -2.969974e-03 -2.586799e-03 -2.769938e-03
## delta.4. -2.718293e-03 -2.298821e-03 -2.498395e-03
## delta.5. -4.051228e-03 -3.821873e-03 -3.933484e-03
## delta.6. -3.588951e-03 -3.197713e-03 -3.382634e-03
## delta.7. -3.879030e-03 -3.684672e-03 -3.779524e-03
## delta.8. -4.119015e-03 -3.975270e-03 -4.046257e-03
## delta.9. -3.750831e-03 -3.637632e-03 -3.693651e-03
## delta.10. -4.088548e-03 -3.894629e-03 -3.989673e-03
## delta.11. -3.756364e-03 -3.595060e-03 -3.674093e-03
## delta.12. -3.611611e-03 -3.482366e-03 -3.545865e-03
## delta.13. -4.320107e-03 -4.105318e-03 -4.209816e-03
## delta.14. -4.097897e-03 -3.946837e-03 -4.021145e-03
## delta.15. -3.980550e-03 -3.808793e-03 -3.892875e-03
## delta.16. -4.056067e-03 -3.886810e-03 -3.970244e-03
## delta.17. -3.837891e-03 -3.730789e-03 -3.783815e-03
## delta.18. -4.166437e-03 -4.030778e-03 -4.097490e-03
## delta.19. -3.668010e-03 -3.564934e-03 -3.615850e-03
## delta.20. -3.587042e-03 -3.497329e-03 -3.541414e-03
## delta.21. -3.400624e-03 -3.292082e-03 -3.346049e-03
## delta.22. -3.974423e-03 -3.823658e-03 -3.897974e-03
## delta.23. -3.749702e-03 -3.652031e-03 -3.700433e-03
## delta.24. -3.351117e-03 -3.254202e-03 -3.302190e-03
## delta.25. -4.264045e-03 -4.100337e-03 -4.180989e-03
## delta.26. -4.122069e-03 -3.945685e-03 -4.032168e-03
## delta.27. -3.631162e-03 -3.536932e-03 -3.583492e-03
## delta.28. -3.966688e-03 -3.860937e-03 -3.913238e-03
## delta.29. -4.278376e-03 -4.086159e-03 -4.179781e-03
## delta.30. -4.747881e-03 -4.358148e-03 -4.543484e-03
## delta.31. -4.150986e-03 -4.012304e-03 -4.080392e-03
## delta.32. -4.613369e-03 -4.136416e-03 -4.360774e-03
## delta.33. -4.320027e-03 -4.140211e-03 -4.228669e-03
## delta.34. -3.480584e-03 -3.371023e-03 -3.425328e-03
## delta.35. -3.604060e-03 -3.491357e-03 -3.546880e-03
## delta.36. -4.438547e-03 -4.213139e-03 -4.321765e-03
## delta.37. -4.011126e-03 -3.890165e-03 -3.949733e-03
## delta.38. -4.046383e-03 -3.898156e-03 -3.970756e-03
## delta.39. -4.120865e-03 -3.969203e-03 -4.043773e-03
## delta.40. -3.825775e-03 -3.611044e-03 -3.715930e-03
## delta.41. -3.908333e-03 -3.764577e-03 -3.835530e-03
## delta.42. -3.463873e-03 -3.365678e-03 -3.414195e-03
## delta.43. -4.320507e-03 -4.100282e-03 -4.207425e-03
## delta.44. -4.819598e-03 -4.399923e-03 -4.598882e-03
## delta.45. -4.029310e-03 -3.903263e-03 -3.965421e-03
## delta.46. -3.567763e-03 -3.456856e-03 -3.511583e-03
## delta.47. -3.645227e-03 -3.541671e-03 -3.592990e-03
## delta.48. -4.309039e-03 -4.080336e-03 -4.191130e-03
## delta.49. -3.715558e-03 -3.564806e-03 -3.638507e-03
## delta.50. -3.663092e-03 -3.465155e-03 -3.561759e-03
## delta.51. -4.471091e-03 -4.157674e-03 -4.307555e-03
## delta.52. -4.593938e-03 -4.240701e-03 -4.410293e-03
## delta.53. -4.154460e-03 -3.995037e-03 -4.072562e-03
## delta.54. -4.328963e-03 -4.090126e-03 -4.206603e-03
## delta.55. -4.244785e-03 -4.029039e-03 -4.134586e-03
## delta.56. -3.888913e-03 -3.554766e-03 -3.714995e-03
## delta.57. -4.794647e-03 -4.448463e-03 -4.614012e-03
## delta.58. -4.593432e-03 -4.245760e-03 -4.413045e-03
## p.1.1. 5.749939e-04 9.602428e-04 7.554488e-04
## p.2.1. 3.131995e-03 5.494800e-03 4.227676e-03
## p.3.1. 2.548631e-03 5.489217e-03 3.877168e-03
## p.4.1. 4.219669e-03 9.750573e-03 6.699385e-03
## p.5.1. 2.906884e-04 4.608558e-04 3.708225e-04
## p.6.1. 7.358186e-04 1.614035e-03 1.136081e-03
## p.7.1. 4.108583e-04 6.071116e-04 5.042648e-04
## p.8.1. 2.536708e-04 3.386202e-04 2.943978e-04
## p.9.1. 5.315555e-04 6.672752e-04 5.972580e-04
## p.10.1. 2.696861e-04 3.981812e-04 3.305846e-04
## p.11.1. 5.256803e-04 7.268444e-04 6.222627e-04
## p.12.1. 7.030785e-04 9.114562e-04 8.040822e-04
## p.13.1. 1.693430e-04 2.607494e-04 2.126423e-04
## p.14.1. 2.646669e-04 3.585291e-04 3.097355e-04
## p.15.1. 3.350469e-04 4.731263e-04 4.011305e-04
## p.16.1. 2.878756e-04 4.044865e-04 3.433503e-04
## p.17.1. 4.462611e-04 5.533940e-04 4.982303e-04
## p.18.1. 2.306134e-04 3.028821e-04 2.655262e-04
## p.19.1. 6.277760e-04 7.721824e-04 6.980831e-04
## p.20.1. 7.386444e-04 8.844757e-04 8.103880e-04
## p.21.1. 1.074039e-03 1.335537e-03 1.200276e-03
## p.22.1. 3.391972e-04 4.592063e-04 3.966882e-04
## p.23.1. 5.327634e-04 6.482515e-04 5.889284e-04
## p.24.1. 1.186280e-03 1.441042e-03 1.310311e-03
## p.25.1. 1.895382e-04 2.633725e-04 2.247577e-04
## p.26.1. 2.521184e-04 3.593598e-04 3.032732e-04
## p.27.1. 6.760041e-04 8.168536e-04 7.448089e-04
## p.28.1. 3.445098e-04 4.260690e-04 3.841264e-04
## p.29.1. 1.841574e-04 2.709842e-04 2.256018e-04
## p.30.1. 7.167902e-05 1.568791e-04 1.102606e-04
## p.31.1. 2.378859e-04 3.143366e-04 2.748372e-04
## p.32.1. 9.392947e-05 2.449535e-04 1.606904e-04
## p.33.1. 1.693701e-04 2.430927e-04 2.043642e-04
## p.34.1. 9.147224e-04 1.139806e-03 1.023656e-03
## p.35.1. 7.138241e-04 8.951459e-04 8.020259e-04
## p.36.1. 1.334728e-04 2.099544e-04 1.699001e-04
## p.37.1. 3.150820e-04 4.017686e-04 3.571354e-04
## p.38.1. 2.935322e-04 3.953701e-04 3.427030e-04
## p.39.1. 2.527294e-04 3.427737e-04 2.959657e-04
## p.40.1. 4.572566e-04 7.038798e-04 5.736094e-04
## p.41.1. 3.873673e-04 5.170781e-04 4.495990e-04
## p.42.1. 9.459398e-04 1.152103e-03 1.046489e-03
## p.43.1. 1.692068e-04 2.634018e-04 2.137256e-04
## p.44.1. 6.205724e-05 1.442461e-04 9.893437e-05
## p.45.1. 3.037772e-04 3.913338e-04 3.461138e-04
## p.46.1. 7.678073e-04 9.593636e-04 8.608960e-04
## p.47.1. 6.571725e-04 8.091164e-04 7.309015e-04
## p.48.1. 1.731517e-04 2.741736e-04 2.209596e-04
## p.49.1. 5.705887e-04 7.723806e-04 6.680907e-04
## p.50.1. 6.340084e-04 9.435087e-04 7.811158e-04
## p.51.1. 1.250224e-04 2.347099e-04 1.759274e-04
## p.52.1. 9.767028e-05 1.986417e-04 1.435782e-04
## p.53.1. 2.362311e-04 3.254346e-04 2.794207e-04
## p.54.1. 1.663558e-04 2.688327e-04 2.143053e-04
## p.55.1. 1.970181e-04 3.039426e-04 2.473618e-04
## p.56.1. 4.027812e-04 7.881140e-04 5.795625e-04
## p.57.1. 6.524867e-05 1.308392e-04 9.527735e-05
## p.58.1. 9.776976e-05 1.966325e-04 1.427305e-04
## p.1.2. 5.728648e-04 9.569328e-04 7.527575e-04
## p.2.2. 3.123056e-03 5.480685e-03 4.216274e-03
## p.3.2. 2.541093e-03 5.475114e-03 3.866558e-03
## p.4.2. 4.208263e-03 9.728402e-03 6.682934e-03
## p.5.2. 2.895134e-04 4.590987e-04 3.693698e-04
## p.6.2. 7.331844e-04 1.608890e-03 1.132271e-03
## p.7.2. 4.092683e-04 6.048801e-04 5.023660e-04
## p.8.2. 2.526283e-04 3.372772e-04 2.932102e-04
## p.9.2. 5.295666e-04 6.648539e-04 5.950583e-04
## p.10.2. 2.685860e-04 3.966341e-04 3.292704e-04
## p.11.2. 5.237104e-04 7.242380e-04 6.199842e-04
## p.12.2. 7.005456e-04 9.082905e-04 8.012402e-04
## p.13.2. 1.686131e-04 2.596814e-04 2.117505e-04
## p.14.2. 2.635849e-04 3.571173e-04 3.084938e-04
## p.15.2. 3.337164e-04 4.713286e-04 3.995741e-04
## p.16.2. 2.867107e-04 4.029180e-04 3.419915e-04
## p.17.2. 4.445524e-04 5.513344e-04 4.963503e-04
## p.18.2. 2.296547e-04 3.016641e-04 2.644414e-04
## p.19.2. 6.254789e-04 7.694366e-04 6.955662e-04
## p.20.2. 7.360015e-04 8.813905e-04 8.075264e-04
## p.21.2. 1.070397e-03 1.331154e-03 1.196273e-03
## p.22.2. 3.378522e-04 4.574546e-04 3.951467e-04
## p.23.2. 5.307705e-04 6.458899e-04 5.867551e-04
## p.24.2. 1.182316e-03 1.436367e-03 1.305999e-03
## p.25.2. 1.887319e-04 2.622950e-04 2.238210e-04
## p.26.2. 2.510815e-04 3.579452e-04 3.020545e-04
## p.27.2. 6.735556e-04 8.139719e-04 7.421475e-04
## p.28.2. 3.431464e-04 4.244278e-04 3.826273e-04
## p.29.2. 1.833713e-04 2.698795e-04 2.246621e-04
## p.30.2. 7.133952e-05 1.561970e-04 1.097630e-04
## p.31.2. 2.369008e-04 3.130783e-04 2.737190e-04
## p.32.2. 9.349718e-05 2.439426e-04 1.599960e-04
## p.33.2. 1.686401e-04 2.420886e-04 2.035028e-04
## p.34.2. 9.115470e-04 1.135974e-03 1.020161e-03
## p.35.2. 7.112579e-04 8.920289e-04 7.991898e-04
## p.36.2. 1.328817e-04 2.090719e-04 1.691687e-04
## p.37.2. 3.138211e-04 4.002093e-04 3.557287e-04
## p.38.2. 2.923472e-04 3.938325e-04 3.413463e-04
## p.39.2. 2.516904e-04 3.414163e-04 2.947725e-04
## p.40.2. 4.555114e-04 7.013444e-04 5.714865e-04
## p.41.2. 3.858569e-04 5.151362e-04 4.478798e-04
## p.42.2. 9.426719e-04 1.148236e-03 1.042927e-03
## p.43.2. 1.684774e-04 2.623243e-04 2.128298e-04
## p.44.2. 6.175889e-05 1.436129e-04 9.848270e-05
## p.45.2. 3.025560e-04 3.898099e-04 3.447452e-04
## p.46.2. 7.650749e-04 9.560561e-04 8.578822e-04
## p.47.2. 6.547829e-04 8.062581e-04 7.282830e-04
## p.48.2. 1.724074e-04 2.730575e-04 2.200372e-04
## p.49.2. 5.684738e-04 7.696343e-04 6.656679e-04
## p.50.2. 6.316917e-04 9.402480e-04 7.783447e-04
## p.51.2. 1.244647e-04 2.337363e-04 1.751736e-04
## p.52.2. 9.722266e-05 1.978012e-04 1.429488e-04
## p.53.2. 2.352520e-04 3.241375e-04 2.782863e-04
## p.54.2. 1.656373e-04 2.677357e-04 2.134074e-04
## p.55.2. 1.961838e-04 3.027209e-04 2.463430e-04
## p.56.2. 4.012185e-04 7.853196e-04 5.774230e-04
## p.57.2. 6.493660e-05 1.302585e-04 9.484026e-05
## p.58.2. 9.732173e-05 1.957996e-04 1.421043e-04
## p.1.3. 5.707437e-04 9.536341e-04 7.500757e-04
## p.2.3. 3.114142e-03 5.466606e-03 4.204903e-03
## p.3.3. 2.533576e-03 5.461046e-03 3.855977e-03
## p.4.3. 4.196887e-03 9.706281e-03 6.666523e-03
## p.5.3. 2.883433e-04 4.573482e-04 3.679228e-04
## p.6.3. 7.305597e-04 1.603762e-03 1.128474e-03
## p.7.3. 4.076845e-04 6.026567e-04 5.004743e-04
## p.8.3. 2.515901e-04 3.359396e-04 2.920273e-04
## p.9.3. 5.275850e-04 6.624414e-04 5.928667e-04
## p.10.3. 2.674904e-04 3.950929e-04 3.279613e-04
## p.11.3. 5.217478e-04 7.216408e-04 6.177140e-04
## p.12.3. 6.980218e-04 9.051359e-04 7.984081e-04
## p.13.3. 1.678864e-04 2.586178e-04 2.108624e-04
## p.14.3. 2.625072e-04 3.557111e-04 3.072571e-04
## p.15.3. 3.323911e-04 4.695376e-04 3.980238e-04
## p.16.3. 2.855504e-04 4.013556e-04 3.406382e-04
## p.17.3. 4.428503e-04 5.492824e-04 4.944775e-04
## p.18.3. 2.287001e-04 3.004510e-04 2.633610e-04
## p.19.3. 6.231903e-04 7.667006e-04 6.930584e-04
## p.20.3. 7.333681e-04 8.783161e-04 8.046748e-04
## p.21.3. 1.066767e-03 1.326784e-03 1.192284e-03
## p.22.3. 3.365125e-04 4.557096e-04 3.936112e-04
## p.23.3. 5.287850e-04 6.435369e-04 5.845899e-04
## p.24.3. 1.178365e-03 1.431707e-03 1.301700e-03
## p.25.3. 1.879290e-04 2.612220e-04 2.228881e-04
## p.26.3. 2.500490e-04 3.565361e-04 3.008406e-04
## p.27.3. 6.711158e-04 8.111003e-04 7.394956e-04
## p.28.3. 3.417884e-04 4.227930e-04 3.811340e-04
## p.29.3. 1.825886e-04 2.687793e-04 2.237263e-04
## p.30.3. 7.100164e-05 1.555179e-04 1.092676e-04
## p.31.3. 2.359196e-04 3.118250e-04 2.726054e-04
## p.32.3. 9.306688e-05 2.429359e-04 1.593045e-04
## p.33.3. 1.679133e-04 2.410886e-04 2.026451e-04
## p.34.3. 9.083827e-04 1.132156e-03 1.016678e-03
## p.35.3. 7.087009e-04 8.889227e-04 7.963638e-04
## p.36.3. 1.322933e-04 2.081931e-04 1.684404e-04
## p.37.3. 3.125652e-04 3.986560e-04 3.543277e-04
## p.38.3. 2.911670e-04 3.923009e-04 3.399951e-04
## p.39.3. 2.506556e-04 3.400643e-04 2.935842e-04
## p.40.3. 4.537728e-04 6.988182e-04 5.693715e-04
## p.41.3. 3.843524e-04 5.132016e-04 4.461672e-04
## p.42.3. 9.394154e-04 1.144383e-03 1.039378e-03
## p.43.3. 1.677512e-04 2.612512e-04 2.119377e-04
## p.44.3. 6.146197e-05 1.429825e-04 9.803308e-05
## p.45.3. 3.013397e-04 3.882920e-04 3.433821e-04
## p.46.3. 7.623523e-04 9.527600e-04 8.548789e-04
## p.47.3. 6.524019e-04 8.034100e-04 7.256739e-04
## p.48.3. 1.716662e-04 2.719459e-04 2.191186e-04
## p.49.3. 5.663667e-04 7.668976e-04 6.632538e-04
## p.50.3. 6.293834e-04 9.369986e-04 7.755835e-04
## p.51.3. 1.239095e-04 2.327668e-04 1.744230e-04
## p.52.3. 9.677709e-05 1.969644e-04 1.423221e-04
## p.53.3. 2.342769e-04 3.228456e-04 2.771565e-04
## p.54.3. 1.649220e-04 2.666432e-04 2.125134e-04
## p.55.3. 1.953529e-04 3.015040e-04 2.453283e-04
## p.56.3. 3.996618e-04 7.825351e-04 5.752914e-04
## p.57.3. 6.462602e-05 1.296804e-04 9.440517e-05
## p.58.3. 9.687576e-05 1.949702e-04 1.414809e-04
## p.1.4. 5.686303e-04 9.503467e-04 7.474036e-04
## p.2.4. 3.105254e-03 5.452563e-03 4.193562e-03
## p.3.4. 2.526081e-03 5.447015e-03 3.845425e-03
## p.4.4. 4.185542e-03 9.684209e-03 6.650153e-03
## p.5.4. 2.871778e-04 4.556044e-04 3.664815e-04
## p.6.4. 7.279443e-04 1.598650e-03 1.124690e-03
## p.7.4. 4.061068e-04 6.004415e-04 4.985897e-04
## p.8.4. 2.505562e-04 3.346072e-04 2.908492e-04
## p.9.4. 5.256109e-04 6.600377e-04 5.906832e-04
## p.10.4. 2.663993e-04 3.935578e-04 3.266575e-04
## p.11.4. 5.197927e-04 7.190530e-04 6.154521e-04
## p.12.4. 6.955071e-04 9.019922e-04 7.955861e-04
## p.13.4. 1.671628e-04 2.575585e-04 2.099780e-04
## p.14.4. 2.614340e-04 3.543105e-04 3.060253e-04
## p.15.4. 3.310710e-04 4.677535e-04 3.964795e-04
## p.16.4. 2.843949e-04 3.997992e-04 3.392902e-04
## p.17.4. 4.411547e-04 5.472381e-04 4.926116e-04
## p.18.4. 2.277494e-04 2.992427e-04 2.622850e-04
## p.19.4. 6.209100e-04 7.639744e-04 6.905596e-04
## p.20.4. 7.307441e-04 8.752524e-04 8.018333e-04
## p.21.4. 1.063149e-03 1.322429e-03 1.188308e-03
## p.22.4. 3.351782e-04 4.539713e-04 3.920817e-04
## p.23.4. 5.268070e-04 6.411925e-04 5.824326e-04
## p.24.4. 1.174427e-03 1.427062e-03 1.297416e-03
## p.25.4. 1.871295e-04 2.601534e-04 2.219592e-04
## p.26.4. 2.490206e-04 3.551326e-04 2.996316e-04
## p.27.4. 6.686850e-04 8.082389e-04 7.368532e-04
## p.28.4. 3.404358e-04 4.211644e-04 3.796466e-04
## p.29.4. 1.818092e-04 2.676835e-04 2.227944e-04
## p.30.4. 7.066536e-05 1.548417e-04 1.087745e-04
## p.31.4. 2.349426e-04 3.105768e-04 2.714963e-04
## p.32.4. 9.263855e-05 2.419334e-04 1.586160e-04
## p.33.4. 1.671896e-04 2.400928e-04 2.017910e-04
## p.34.4. 9.052294e-04 1.128350e-03 1.013206e-03
## p.35.4. 7.061531e-04 8.858273e-04 7.935478e-04
## p.36.4. 1.317075e-04 2.073180e-04 1.677152e-04
## p.37.4. 3.113144e-04 3.971088e-04 3.529321e-04
## p.38.4. 2.899915e-04 3.907752e-04 3.386492e-04
## p.39.4. 2.496250e-04 3.387177e-04 2.924006e-04
## p.40.4. 4.520409e-04 6.963010e-04 5.672643e-04
## p.41.4. 3.828537e-04 5.112742e-04 4.444612e-04
## p.42.4. 9.361700e-04 1.140542e-03 1.035840e-03
## p.43.4. 1.670281e-04 2.601825e-04 2.110494e-04
## p.44.4. 6.116648e-05 1.423549e-04 9.758552e-05
## p.45.4. 3.001283e-04 3.867799e-04 3.420243e-04
## p.46.4. 7.596393e-04 9.494752e-04 8.518861e-04
## p.47.4. 6.500296e-04 8.005719e-04 7.230742e-04
## p.48.4. 1.709282e-04 2.708388e-04 2.182039e-04
## p.49.4. 5.642674e-04 7.641708e-04 6.608485e-04
## p.50.4. 6.270836e-04 9.337604e-04 7.728321e-04
## p.51.4. 1.233568e-04 2.318012e-04 1.736756e-04
## p.52.4. 9.633356e-05 1.961310e-04 1.416982e-04
## p.53.4. 2.333058e-04 3.215588e-04 2.760313e-04
## p.54.4. 1.642097e-04 2.655551e-04 2.116230e-04
## p.55.4. 1.945256e-04 3.002921e-04 2.443178e-04
## p.56.4. 3.981112e-04 7.797605e-04 5.731676e-04
## p.57.4. 6.431692e-05 1.291049e-04 9.397207e-05
## p.58.4. 9.643183e-05 1.941443e-04 1.408602e-04
## p.1.5. 5.665248e-04 9.470707e-04 7.447409e-04
## p.2.5. 3.096391e-03 5.438556e-03 4.182252e-03
## p.3.5. 2.518609e-03 5.433019e-03 3.834902e-03
## p.4.5. 4.174227e-03 9.662187e-03 6.633822e-03
## p.5.5. 2.860171e-04 4.538672e-04 3.650458e-04
## p.6.5. 7.253384e-04 1.593554e-03 1.120919e-03
## p.7.5. 4.045352e-04 5.982345e-04 4.967122e-04
## p.8.5. 2.495265e-04 3.332802e-04 2.896758e-04
## p.9.5. 5.236441e-04 6.576426e-04 5.885077e-04
## p.10.5. 2.653126e-04 3.920286e-04 3.253589e-04
## p.11.5. 5.178448e-04 7.164745e-04 6.131985e-04
## p.12.5. 6.930015e-04 8.988594e-04 7.927740e-04
## p.13.5. 1.664423e-04 2.565036e-04 2.090973e-04
## p.14.5. 2.603651e-04 3.529153e-04 3.047985e-04
## p.15.5. 3.297563e-04 4.659761e-04 3.949412e-04
## p.16.5. 2.832440e-04 3.982489e-04 3.379475e-04
## p.17.5. 4.394656e-04 5.452014e-04 4.907529e-04
## p.18.5. 2.268027e-04 2.980393e-04 2.612134e-04
## p.19.5. 6.186381e-04 7.612578e-04 6.880698e-04
## p.20.5. 7.281295e-04 8.721994e-04 7.990018e-04
## p.21.5. 1.059544e-03 1.318089e-03 1.184345e-03
## p.22.5. 3.338492e-04 4.522395e-04 3.905581e-04
## p.23.5. 5.248364e-04 6.388566e-04 5.802833e-04
## p.24.5. 1.170503e-03 1.422432e-03 1.293146e-03
## p.25.5. 1.863334e-04 2.590891e-04 2.210341e-04
## p.26.5. 2.479965e-04 3.537346e-04 2.984274e-04
## p.27.5. 6.662629e-04 8.053876e-04 7.342202e-04
## p.28.5. 3.390885e-04 4.195422e-04 3.781650e-04
## p.29.5. 1.810332e-04 2.665922e-04 2.218664e-04
## p.30.5. 7.033066e-05 1.541684e-04 1.082836e-04
## p.31.5. 2.339696e-04 3.093335e-04 2.703917e-04
## p.32.5. 9.221220e-05 2.409349e-04 1.579305e-04
## p.33.5. 1.664690e-04 2.391010e-04 2.009406e-04
## p.34.5. 9.020869e-04 1.124557e-03 1.009747e-03
## p.35.5. 7.036145e-04 8.827427e-04 7.907417e-04
## p.36.5. 1.311243e-04 2.064465e-04 1.669932e-04
## p.37.5. 3.100685e-04 3.955676e-04 3.515420e-04
## p.38.5. 2.888208e-04 3.892555e-04 3.373086e-04
## p.39.5. 2.485987e-04 3.373763e-04 2.912218e-04
## p.40.5. 4.503156e-04 6.937929e-04 5.651649e-04
## p.41.5. 3.813609e-04 5.093541e-04 4.427617e-04
## p.42.5. 9.329359e-04 1.136714e-03 1.032314e-03
## p.43.5. 1.663081e-04 2.591181e-04 2.101648e-04
## p.44.5. 6.087241e-05 1.417300e-04 9.714000e-05
## p.45.5. 2.989218e-04 3.852737e-04 3.406719e-04
## p.46.5. 7.569359e-04 9.462018e-04 8.489037e-04
## p.47.5. 6.476660e-04 7.977438e-04 7.204837e-04
## p.48.5. 1.701933e-04 2.697363e-04 2.172930e-04
## p.49.5. 5.621759e-04 7.614536e-04 6.584519e-04
## p.50.5. 6.247921e-04 9.305334e-04 7.700904e-04
## p.51.5. 1.228066e-04 2.308397e-04 1.729315e-04
## p.52.5. 9.589207e-05 1.953012e-04 1.410770e-04
## p.53.5. 2.323388e-04 3.202771e-04 2.749106e-04
## p.54.5. 1.635005e-04 2.644714e-04 2.107364e-04
## p.55.5. 1.937018e-04 2.990850e-04 2.433114e-04
## p.56.5. 3.965666e-04 7.769957e-04 5.710517e-04
## p.57.5. 6.400930e-05 1.285319e-04 9.354097e-05
## p.58.5. 9.598994e-05 1.933219e-04 1.402422e-04
## p.1.6. 5.644271e-04 9.438060e-04 7.420878e-04
## p.2.6. 3.087553e-03 5.424585e-03 4.170972e-03
## p.3.6. 2.511159e-03 5.419059e-03 3.824407e-03
## p.4.6. 4.162943e-03 9.640215e-03 6.617532e-03
## p.5.6. 2.848610e-04 4.521367e-04 3.636157e-04
## p.6.6. 7.227417e-04 1.588474e-03 1.117160e-03
## p.7.6. 4.029696e-04 5.960356e-04 4.948418e-04
## p.8.6. 2.485011e-04 3.319584e-04 2.885072e-04
## p.9.6. 5.216847e-04 6.552563e-04 5.863403e-04
## p.10.6. 2.642304e-04 3.905054e-04 3.240654e-04
## p.11.6. 5.159042e-04 7.139052e-04 6.109531e-04
## p.12.6. 6.905049e-04 8.957375e-04 7.899718e-04
## p.13.6. 1.657249e-04 2.554530e-04 2.082204e-04
## p.14.6. 2.593006e-04 3.515256e-04 3.035766e-04
## p.15.6. 3.284467e-04 4.642055e-04 3.934088e-04
## p.16.6. 2.820978e-04 3.967046e-04 3.366102e-04
## p.17.6. 4.377829e-04 5.431723e-04 4.889011e-04
## p.18.6. 2.258599e-04 2.968408e-04 2.601462e-04
## p.19.6. 6.163745e-04 7.585508e-04 6.855890e-04
## p.20.6. 7.255242e-04 8.691570e-04 7.961803e-04
## p.21.6. 1.055950e-03 1.313762e-03 1.180395e-03
## p.22.6. 3.325254e-04 4.505144e-04 3.890404e-04
## p.23.6. 5.228731e-04 6.365292e-04 5.781420e-04
## p.24.6. 1.166592e-03 1.417818e-03 1.288890e-03
## p.25.6. 1.855407e-04 2.580292e-04 2.201129e-04
## p.26.6. 2.469766e-04 3.523422e-04 2.972281e-04
## p.27.6. 6.638495e-04 8.025463e-04 7.315966e-04
## p.28.6. 3.377466e-04 4.179261e-04 3.766891e-04
## p.29.6. 1.802604e-04 2.655054e-04 2.209423e-04
## p.30.6. 6.999756e-05 1.534981e-04 1.077949e-04
## p.31.6. 2.330006e-04 3.080953e-04 2.692917e-04
## p.32.6. 9.178781e-05 2.399406e-04 1.572480e-04
## p.33.6. 1.657515e-04 2.381134e-04 2.000937e-04
## p.34.6. 8.989554e-04 1.120777e-03 1.006299e-03
## p.35.6. 7.010849e-04 8.796688e-04 7.879455e-04
## p.36.6. 1.305436e-04 2.055787e-04 1.662743e-04
## p.37.6. 3.088277e-04 3.940324e-04 3.501574e-04
## p.38.6. 2.876548e-04 3.877416e-04 3.359733e-04
## p.39.6. 2.475766e-04 3.360403e-04 2.900478e-04
## p.40.6. 4.485968e-04 6.912939e-04 5.630733e-04
## p.41.6. 3.798738e-04 5.074411e-04 4.410686e-04
## p.42.6. 9.297129e-04 1.132899e-03 1.028801e-03
## p.43.6. 1.655913e-04 2.580581e-04 2.092839e-04
## p.44.6. 6.057975e-05 1.411078e-04 9.669652e-05
## p.45.6. 2.977201e-04 3.837734e-04 3.393248e-04
## p.46.6. 7.542422e-04 9.429396e-04 8.459319e-04
## p.47.6. 6.453109e-04 7.949257e-04 7.179025e-04
## p.48.6. 1.694617e-04 2.686382e-04 2.163858e-04
## p.49.6. 5.600921e-04 7.587460e-04 6.560640e-04
## p.50.6. 6.225091e-04 9.273175e-04 7.673584e-04
## p.51.6. 1.222588e-04 2.298822e-04 1.721905e-04
## p.52.6. 9.545260e-05 1.944749e-04 1.404585e-04
## p.53.6. 2.313758e-04 3.190005e-04 2.737945e-04
## p.54.6. 1.627943e-04 2.633922e-04 2.098535e-04
## p.55.6. 1.928815e-04 2.978827e-04 2.423092e-04
## p.56.6. 3.950280e-04 7.742407e-04 5.689436e-04
## p.57.6. 6.370315e-05 1.279615e-04 9.311184e-05
## p.58.6. 9.555007e-05 1.925030e-04 1.396270e-04
## p.1.7. 5.623372e-04 9.405525e-04 7.394441e-04
## p.2.7. 3.078741e-03 5.410649e-03 4.159722e-03
## p.3.7. 2.503730e-03 5.405135e-03 3.813941e-03
## p.4.7. 4.151689e-03 9.618292e-03 6.601281e-03
## p.5.7. 2.837097e-04 4.504128e-04 3.621913e-04
## p.6.7. 7.201543e-04 1.583411e-03 1.113414e-03
## p.7.7. 4.014102e-04 5.938447e-04 4.929784e-04
## p.8.7. 2.474799e-04 3.306418e-04 2.873433e-04
## p.9.7. 5.197326e-04 6.528786e-04 5.841808e-04
## p.10.7. 2.631525e-04 3.889881e-04 3.227771e-04
## p.11.7. 5.139709e-04 7.113451e-04 6.087160e-04
## p.12.7. 6.880173e-04 8.926264e-04 7.871796e-04
## p.13.7. 1.650107e-04 2.544067e-04 2.073471e-04
## p.14.7. 2.582405e-04 3.501414e-04 3.023597e-04
## p.15.7. 3.271423e-04 4.624417e-04 3.918824e-04
## p.16.7. 2.809563e-04 3.951663e-04 3.352781e-04
## p.17.7. 4.361067e-04 5.411507e-04 4.870563e-04
## p.18.7. 2.249211e-04 2.956470e-04 2.590833e-04
## p.19.7. 6.141192e-04 7.558535e-04 6.831171e-04
## p.20.7. 7.229283e-04 8.661252e-04 7.933687e-04
## p.21.7. 1.052369e-03 1.309450e-03 1.176459e-03
## p.22.7. 3.312068e-04 4.487958e-04 3.875286e-04
## p.23.7. 5.209172e-04 6.342103e-04 5.760085e-04
## p.24.7. 1.162693e-03 1.413218e-03 1.284648e-03
## p.25.7. 1.847514e-04 2.569737e-04 2.191955e-04
## p.26.7. 2.459609e-04 3.509551e-04 2.960336e-04
## p.27.7. 6.614450e-04 7.997150e-04 7.289824e-04
## p.28.7. 3.364099e-04 4.163163e-04 3.752190e-04
## p.29.7. 1.794910e-04 2.644230e-04 2.200220e-04
## p.30.7. 6.966603e-05 1.528307e-04 1.073084e-04
## p.31.7. 2.320357e-04 3.068619e-04 2.681961e-04
## p.32.7. 9.136537e-05 2.389504e-04 1.565684e-04
## p.33.7. 1.650371e-04 2.371298e-04 1.992503e-04
## p.34.7. 8.958348e-04 1.117009e-03 1.002863e-03
## p.35.7. 6.985645e-04 8.766056e-04 7.851592e-04
## p.36.7. 1.299656e-04 2.047146e-04 1.655584e-04
## p.37.7. 3.075918e-04 3.925031e-04 3.487783e-04
## p.38.7. 2.864935e-04 3.862337e-04 3.346433e-04
## p.39.7. 2.465588e-04 3.347096e-04 2.888784e-04
## p.40.7. 4.468846e-04 6.888038e-04 5.609893e-04
## p.41.7. 3.783926e-04 5.055354e-04 4.393821e-04
## p.42.7. 9.265010e-04 1.129097e-03 1.025299e-03
## p.43.7. 1.648775e-04 2.570024e-04 2.084067e-04
## p.44.7. 6.028850e-05 1.404884e-04 9.625506e-05
## p.45.7. 2.965233e-04 3.822789e-04 3.379831e-04
## p.46.7. 7.515581e-04 9.396887e-04 8.429704e-04
## p.47.7. 6.429644e-04 7.921175e-04 7.153306e-04
## p.48.7. 1.687332e-04 2.675446e-04 2.154825e-04
## p.49.7. 5.580161e-04 7.560481e-04 6.536847e-04
## p.50.7. 6.202344e-04 9.241127e-04 7.646362e-04
## p.51.7. 1.217134e-04 2.289286e-04 1.714527e-04
## p.52.7. 9.501515e-05 1.936521e-04 1.398428e-04
## p.53.7. 2.304168e-04 3.177291e-04 2.726830e-04
## p.54.7. 1.620912e-04 2.623174e-04 2.089743e-04
## p.55.7. 1.920646e-04 2.966853e-04 2.413112e-04
## p.56.7. 3.934953e-04 7.714955e-04 5.668433e-04
## p.57.7. 6.339847e-05 1.273936e-04 9.268468e-05
## p.58.7. 9.511221e-05 1.916876e-04 1.390144e-04
## p.1.8. 5.602550e-04 9.373103e-04 7.368098e-04
## p.2.8. 3.069953e-03 5.396749e-03 4.148503e-03
## p.3.8. 2.496324e-03 5.391246e-03 3.803504e-03
## p.4.8. 4.140466e-03 9.596419e-03 6.585070e-03
## p.5.8. 2.825629e-04 4.486954e-04 3.607724e-04
## p.6.8. 7.175762e-04 1.578364e-03 1.109681e-03
## p.7.8. 3.998567e-04 5.916619e-04 4.911220e-04
## p.8.8. 2.464629e-04 3.293304e-04 2.861841e-04
## p.9.8. 5.177879e-04 6.505095e-04 5.820292e-04
## p.10.8. 2.620791e-04 3.874766e-04 3.214939e-04
## p.11.8. 5.120449e-04 7.087941e-04 6.064870e-04
## p.12.8. 6.855386e-04 8.895261e-04 7.843972e-04
## p.13.8. 1.642994e-04 2.533647e-04 2.064774e-04
## p.14.8. 2.571847e-04 3.487627e-04 3.011475e-04
## p.15.8. 3.258431e-04 4.606845e-04 3.903620e-04
## p.16.8. 2.798193e-04 3.936340e-04 3.339513e-04
## p.17.8. 4.344369e-04 5.391366e-04 4.852185e-04
## p.18.8. 2.239861e-04 2.944581e-04 2.580248e-04
## p.19.8. 6.118721e-04 7.531657e-04 6.806541e-04
## p.20.8. 7.203416e-04 8.631040e-04 7.905671e-04
## p.21.8. 1.048800e-03 1.305152e-03 1.172535e-03
## p.22.8. 3.298935e-04 4.470838e-04 3.860227e-04
## p.23.8. 5.189685e-04 6.318999e-04 5.738829e-04
## p.24.8. 1.158808e-03 1.408633e-03 1.280420e-03
## p.25.8. 1.839655e-04 2.559224e-04 2.182819e-04
## p.26.8. 2.449494e-04 3.495736e-04 2.948439e-04
## p.27.8. 6.590491e-04 7.968937e-04 7.263775e-04
## p.28.8. 3.350786e-04 4.147127e-04 3.737547e-04
## p.29.8. 1.787248e-04 2.633450e-04 2.191055e-04
## p.30.8. 6.933607e-05 1.521662e-04 1.068241e-04
## p.31.8. 2.310747e-04 3.056336e-04 2.671049e-04
## p.32.8. 9.094488e-05 2.379643e-04 1.558917e-04
## p.33.8. 1.643258e-04 2.361503e-04 1.984105e-04
## p.34.8. 8.927250e-04 1.113254e-03 9.994387e-04
## p.35.8. 6.960531e-04 8.735531e-04 7.823827e-04
## p.36.8. 1.293901e-04 2.038541e-04 1.648457e-04
## p.37.8. 3.063609e-04 3.909798e-04 3.474046e-04
## p.38.8. 2.853370e-04 3.847316e-04 3.333185e-04
## p.39.8. 2.455451e-04 3.333841e-04 2.877138e-04
## p.40.8. 4.451790e-04 6.863227e-04 5.589132e-04
## p.41.8. 3.769172e-04 5.036368e-04 4.377020e-04
## p.42.8. 9.233003e-04 1.125307e-03 1.021809e-03
## p.43.8. 1.641668e-04 2.559511e-04 2.075332e-04
## p.44.8. 5.999865e-05 1.398717e-04 9.581562e-05
## p.45.8. 2.953313e-04 3.807903e-04 3.366466e-04
## p.46.8. 7.488835e-04 9.364490e-04 8.400192e-04
## p.47.8. 6.406264e-04 7.893193e-04 7.127679e-04
## p.48.8. 1.680078e-04 2.664554e-04 2.145829e-04
## p.49.8. 5.559478e-04 7.533597e-04 6.513141e-04
## p.50.8. 6.179679e-04 9.209190e-04 7.619236e-04
## p.51.8. 1.211705e-04 2.279790e-04 1.707181e-04
## p.52.8. 9.457969e-05 1.928328e-04 1.392297e-04
## p.53.8. 2.294617e-04 3.164627e-04 2.715759e-04
## p.54.8. 1.613912e-04 2.612469e-04 2.080988e-04
## p.55.8. 1.912512e-04 2.954927e-04 2.403172e-04
## p.56.8. 3.919686e-04 7.687600e-04 5.647507e-04
## p.57.8. 6.309524e-05 1.268282e-04 9.225948e-05
## p.58.8. 9.467636e-05 1.908756e-04 1.384045e-04The 95% confidence interval for \(\beta\) is (-0.004096, 0.0002897). Since \(\beta\) is not considered signifant in the 5% level, we can proceed with Model 2 instead, in which the overall time trend is removed from the model.
Crime Model 1 Visualizations
#Alter the data frame for visualizations
zip_data <- data.frame(zipcode=c(1:58), zip=rownames(Crime_count))
crime_viz_1 <- crime_data_1 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "delta")) %>%
filter(grepl("delta", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))Map of Crime Differential Growth by Zipcode (Using Model 1)
zip <- readOGR(dsn="Boundaries - ZIP Codes", layer="geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a")## OGR data source with driver: ESRI Shapefile
## Source: "Boundaries - ZIP Codes", layer: "geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a"
## with 61 features
## It has 4 fieldszip_df <- fortify(zip, region = "zip")
crime_model_1_delta <- left_join(zip_df, crime_viz_1, by="id")
head(crime_model_1_delta)## long lat order hole piece id group zipcode delta
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 -0.003580096
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 -0.003580096
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 -0.003580096
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 -0.003580096
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 -0.003580096
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 -0.003580096
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601model_1_map <- ggplot()
model_1_map <- model_1_map + geom_polygon(data = crime_model_1_delta, aes(x=long, y=lat, group=group, fill=delta), color = "black", size=0.2)
model_1_map <- model_1_map + coord_map()
model_1_map <- model_1_map + labs(x="", y="", title=paste("Crime Differential Growth by Zipcode"), fill="total")
model_1_map <- model_1_map + cleanup
print(model_1_map)
Model 3 – Use demographics information
Notation
i: zipcodes in Chicago (58 total in our data)
j: years (from 2010 to 2017)
y[i, j]: observed count of violent crimes in Chicago by year and zipcode
p[i, j]: probability of being a victim in location i and time j
n[i, j]: the population in zipcode i in time j.
time[j]: the time in years
\(\delta\)[i]: the spatial temporal interaction term
\(\sigma\)[i]: the income trend
income[i]: income for each zipcode in Chicago
Specify the Model
This model resembles Model 2 but now our prediction of the crime rate accounts for the different income levels in various Chicago zipcodes.
Data
\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \delta_{i} * time_{j} + \sigma_{i} * income_{i}\]
Priors
\[\delta_{i} \sim N(0, prec.delta)\] \[\sigma_{i} \sim N(0, prec.sigma)\]
Hyperprior
\[prec.u \sim Gamma(0.5, 0.0005)\] \[prec.delta \sim Gamma(0.5, 0.0005)\] \[prec.sigma \sim Gamma(0.5, 0.0005)\]
#Specify the model
crime_model_3 <- "model{
#Data
for(i in 1:58) {
for(j in 1:8){
y[i,j] ~ dbin(p[i,j], n[i])
logit(p[i,j]) = delta[i] * time[j] + sigma[i] * income[i]
}
}
#Priors
for(i in 1:58){
delta[i] ~ dnorm(0, prec.delta)
sigma[i] ~ dnorm(0, prec.sigma)
}
#Hyperpriors
prec.delta ~ dgamma(0.5, 0.0005)
prec.sigma ~ dgamma(0.5, 0.0005)
}"
#set up an algorithm to simulate the posterior by
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_3 <- jags.model(textConnection(crime_model_3), data=list(y=Crime_count, n=Population.vector, time=time, income=Model_Data$per_capita_income),
inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 464
## Unobserved stochastic nodes: 118
## Total graph size: 2564
##
## Initializing model#simulate a sample from the posterior
#note that we specify both mu and tau variables
crime_sim_3 <- coda.samples(crime_jags_3, variable.names=c("p", "delta", "sigma"), n.iter=40000)
#store the samples in a data frame:
crime_data_3 <- data.frame(crime_sim_3[[1]])Crime Model 3 Visualizations
crime_viz_3_delta <- crime_data_3 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "delta")) %>%
filter(grepl("delta", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
crime_model_3_delta <- left_join(zip_df, crime_viz_3_delta, by="id")
head(crime_model_3_delta)## long lat order hole piece id group zipcode delta
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 -0.1297566
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 -0.1297566
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 -0.1297566
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 -0.1297566
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 -0.1297566
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 -0.1297566
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601crime_viz_3_sigma <- crime_data_3 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "sigma")) %>%
filter(grepl("sigma", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
crime_model_3_sigma <- left_join(zip_df, crime_viz_3_sigma, by="id")
head(crime_model_3_sigma)## long lat order hole piece id group zipcode sigma
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 0.002849706
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 0.002849706
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 0.002849706
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 0.002849706
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 0.002849706
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 0.002849706
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601Map of Crime Differential Growth by Zipcode (Using Model 3)
model_3_map <- ggplot()
model_3_map <- model_3_map + geom_polygon(data = crime_model_3_delta, aes(x=long, y=lat, group=group, fill=delta), color = "black", size=0.2)
model_3_map <- model_3_map + coord_map()
model_3_map <- model_3_map + labs(x="", y="", title=paste("Crime Differential Growth by Zipcode"), fill="total")
model_3_map <- model_3_map + cleanup
print(model_3_map)
Income Differential Growth by Zipcode (Using Model 3)
model_3_map <- ggplot()
model_3_map <- model_3_map + geom_polygon(data = crime_model_3_sigma, aes(x=long, y=lat, group=group, fill=sigma), color = "black", size=0.2)
model_3_map <- model_3_map + coord_map()
model_3_map <- model_3_map + labs(x="", y="", title=paste("Income Differential Growth by Zipcode"), fill="total")
model_3_map <- model_3_map + cleanup
print(model_3_map)
In locations where the income differential growth increases, the crime differential growth seems to decrease.
Running mean plot of alpha
The running mean of \(\alpha\) appears to stabilize after 40,000 iterations.
running_mean_plot(crime_data_3$delta.1.)
running_mean_plot(crime_data_3$delta.3.)
running_mean_plot(crime_data_3$sigma.1.)
running_mean_plot(crime_data_3$sigma.3.)
Model 4 – Use demographics information with race
Notation
i: zipcodes in Chicago (58 total in our data)
j: years (from 2010 to 2017)
y[i, j]: observed count of violent crimes in Chicago by year and zipcode
p[i, j]: probability of being a victim in location i and time j
n[i, j]: the population in zipcode i in time j.
\(\alpha\): the mean crime rate
u[i]: the unstructured spatial random effect for each zipcode in Chicago.
s[i]: the structured spatial random effect for each zipcode in Chicago.
time[j]: the time in years
\(\delta\)[i]: the spatial temporal interaction term
\(\sigma\)[i]: the income trend
income[i]: income for each zipcode in Chicago
e[i]: the weight in the relationship between crime rate and percentage of nonewhite people in Chicago
percent_nonwhite[i]: the percentage of nonwhite people in zipcode i.
Specify the Model
This model resembles Model 3 but now our prediction of the crime rate accounts for both income levels and the percentage of nonwhites in the area.
Data
\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \alpha + u_{i} + \delta_{i} * time_{j} + \sigma_{i} * income_{i} + \rho_{i} * percentNonwhite_{i}\]
Priors
\[\alpha \sim N(0, 1000^2)\] \[\beta \sim N(0, 1000^2)\] \[u_{i} \sim N(0, prec.u)\] \[\delta_{i} \sim N(0, prec.delta)\] \[\sigma_{i} \sim N(0, prec.sigma)\] \[\rho_{i} \sim N(0, 0.01)\]
Hyperprior
\[prec.u \sim Gamma(0.5, 0.0005)\] \[prec.delta \sim Gamma(0.5, 0.0005)\] \[prec.sigma \sim Gamma(0.5, 0.0005)\]
#Specify the model
crime_model_4 <- "model{
#Data
for(i in 1:58) {
for(j in 1:8){
y[i,j] ~ dbin(p[i,j], n[i])
logit(p[i,j]) = delta[i] * time[j] + sigma[i]*income[i] + rho[i]*percent_nonwhite[i]
}
}
#Priors
for(i in 1:58){
delta[i] ~ dnorm(0, prec.delta)
sigma[i] ~ dnorm(0, prec.sigma)
rho[i] ~ dnorm(0, 0.01)
}
#Hyperpriors
prec.delta ~ dgamma(0.5, 0.0005)
prec.sigma ~ dgamma(0.5, 0.0005)
}"
#set up an algorithm to simulate the posterior by
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_4 <- jags.model(textConnection(crime_model_4), data=list(y=Crime_count, n=Population.vector, time=time, income=Model_Data$per_capita_income, percent_nonwhite=100 - Model_Data$percent_white),
inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 464
## Unobserved stochastic nodes: 176
## Total graph size: 3144
##
## Initializing model#simulate a sample from the posterior
#note that we specify both mu and tau variables
crime_sim_4 <- coda.samples(crime_jags_4, variable.names=c("p", "delta", "sigma", "rho"), n.iter=40000)
#store the samples in a data frame:
crime_data_4 <- data.frame(crime_sim_4[[1]])Running Mean Plots
running_mean_plot(crime_data_4$delta.4.)
running_mean_plot(crime_data_4$sigma.4.)
running_mean_plot(crime_data_4$rho.4.)
Crime Model 4 Visualizations
crime_viz_4_delta <- crime_data_4 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "delta")) %>%
filter(grepl("delta", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
crime_model_4_delta <- left_join(zip_df, crime_viz_4_delta, by="id")
head(crime_model_4_delta)## long lat order hole piece id group zipcode delta
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 -0.02008539
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 -0.02008539
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 -0.02008539
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 -0.02008539
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 -0.02008539
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 -0.02008539
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601crime_viz_4_sigma <- crime_data_4 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "sigma")) %>%
filter(grepl("sigma", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
crime_model_4_sigma <- left_join(zip_df, crime_viz_4_sigma, by="id")
head(crime_model_4_sigma)## long lat order hole piece id group zipcode sigma
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 0.0001368507
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 0.0001368507
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 0.0001368507
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 0.0001368507
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 0.0001368507
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 0.0001368507
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601crime_viz_4_rho <- crime_data_4 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "rho")) %>%
filter(grepl("rho", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
crime_model_4_rho <- left_join(zip_df, crime_viz_4_rho, by="id")
head(crime_model_4_rho)## long lat order hole piece id group zipcode rho
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 0.5684385
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 0.5684385
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 0.5684385
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 0.5684385
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 0.5684385
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 0.5684385
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601Map of Crime Differential Growth by Zipcode (Using Model 4)
model_4_map <- ggplot()
model_4_map <- model_4_map + geom_polygon(data = crime_model_4_delta, aes(x=long, y=lat, group=group, fill=delta), color = "black", size=0.2)
model_4_map <- model_4_map + coord_map()
model_4_map <- model_4_map + labs(x="", y="", title=paste("Crime Differential Growth by Zipcode"), fill="total")
model_4_map <- model_4_map + cleanup
print(model_4_map)
Map of Income Differential Growth by Zipcode (Using Model 4)
model_4_map_sigma <- ggplot()
model_4_map_sigma <- model_4_map_sigma + geom_polygon(data = crime_model_4_sigma, aes(x=long, y=lat, group=group, fill=sigma), color = "black", size=0.2)
model_4_map_sigma <- model_4_map_sigma + coord_map()
model_4_map_sigma <- model_4_map_sigma + labs(x="", y="", title=paste("Income Differential Growth by Zipcode"), fill="total")
model_4_map_sigma <- model_4_map_sigma + cleanup
print(model_4_map_sigma)
Map of Percentage Nonwhite Differential Growth by Zipcode (Using Model 4)
model_4_map_rho <- ggplot()
model_4_map_rho <- model_4_map_rho + geom_polygon(data = crime_model_4_rho, aes(x=long, y=lat, group=group, fill=rho), color = "black", size=0.2)
model_4_map_rho <- model_4_map_rho + coord_map()
model_4_map_rho <- model_4_map_rho + labs(x="", y="", title=paste("Percentage Nonwhite Differential Growth by Zipcode"), fill="total")
model_4_map_rho <- model_4_map_rho + cleanup
print(model_4_map_rho)
Police Traffic Stops model
Cleaning the data
Data frame for the total population by race and zipcode
#Number of zipcodes: 62; Number of races: 7
Police <- read.csv("Police_stops_2016_2017_sample.csv")
Police.r <- Police %>%
filter(!(RACE_CODE_CD=="I" | RACE_CODE_CD=="P" | RACE_CODE_CD=="WHI") ) %>%
mutate(RACE_CODE_CD = as.character(RACE_CODE_CD)) %>%
mutate(RACE_CODE_CD=replace(RACE_CODE_CD, RACE_CODE_CD=="WWH", "HISPANIC"))
vars = c("percent_white","percent_black","percent_asian","percent_hispanic")
Population.race <- Model_Data %>%
mutate_at(vars, function(x){trunc(Model_Data$total_population*x/100)}) %>%
select(zip, percent_asian,percent_black,percent_hispanic, percent_white) %>%
unique() %>%
remove_rownames %>%
column_to_rownames(var='zip')
head(Population.race)## percent_asian percent_black percent_hispanic percent_white
## 60601 1751 1070 583 6131
## 60602 321 29 175 907
## 60603 264 35 44 538
## 60604 66 0 12 332
## 60605 4245 4494 1248 14483
## 60606 362 55 251 2091Data frame for the number of arrests by zipcode and race
valid_zips <- unique(rownames(Population.race))
Arrests <- Police.r %>%
filter(ZIP_CD %in% valid_zips) %>%
group_by(ZIP_CD, RACE_CODE_CD) %>%
dplyr::summarise(total.arrest = as.integer(sum(ENFORCEMENT_ACTION_TAKEN_I))) %>%
spread(key=ZIP_CD, value=total.arrest) %>%
as.data.frame() %>%
remove_rownames %>%
column_to_rownames(var='RACE_CODE_CD')
Arrests[is.na(Arrests)] <- as.integer(0)
head(Arrests)## 60601 60602 60603 60604 60605 60606 60607 60608 60609 60610 60611
## API 0 0 0 0 0 0 0 0 1 1 0
## BLK 7 2 2 3 5 2 13 13 87 41 3
## HISPANIC 3 0 1 0 1 0 1 39 43 1 0
## WHT 2 0 0 0 4 0 5 7 8 3 0
## 60612 60613 60614 60615 60616 60617 60618 60619 60620 60621 60622
## API 0 1 0 0 3 0 0 0 1 0 0
## BLK 76 9 7 39 16 100 3 88 103 97 7
## HISPANIC 6 6 3 0 1 19 28 0 2 0 10
## WHT 9 7 4 0 3 3 11 0 3 0 6
## 60623 60624 60625 60626 60628 60629 60630 60631 60632 60633 60634
## API 0 2 1 0 4 0 0 0 0 0 0
## BLK 147 352 6 31 121 40 1 1 5 0 3
## HISPANIC 27 9 21 7 0 27 3 1 40 1 8
## WHT 7 12 4 2 1 4 7 3 5 0 4
## 60636 60637 60638 60639 60640 60641 60642 60643 60644 60645 60646
## API 1 0 0 0 1 0 0 0 0 0 0
## BLK 128 92 1 9 29 4 1 25 211 4 0
## HISPANIC 2 1 1 25 7 13 3 0 3 7 0
## WHT 2 0 4 4 7 8 1 2 9 3 1
## 60647 60649 60651 60652 60653 60654 60655 60656 60657 60659 60660
## API 0 0 0 0 0 0 0 0 0 3 2
## BLK 4 63 96 7 42 6 0 0 10 1 11
## HISPANIC 23 1 29 5 0 0 1 1 4 4 1
## WHT 3 1 6 2 0 1 2 0 2 6 2
## 60661 60707 60827
## API 0 0 0
## BLK 1 2 7
## HISPANIC 0 3 0
## WHT 1 2 0Data frame for the number of police stops by zipcode and race
Stops <- Police.r %>%
filter(ZIP_CD %in% valid_zips) %>%
group_by(ZIP_CD, RACE_CODE_CD) %>%
dplyr::summarise(total.stop = n()) %>%
spread(key=ZIP_CD, value=total.stop) %>%
as.data.frame() %>%
remove_rownames %>%
column_to_rownames(var='RACE_CODE_CD')
Stops[is.na(Stops)] <- 0
head(Stops)## 60601 60602 60603 60604 60605 60606 60607 60608 60609 60610 60611
## API 0 0 0 0 0 0 0 0 2 2 0
## BLK 15 4 6 10 16 6 36 50 291 121 9
## HISPANIC 4 0 1 1 4 0 2 192 187 7 3
## WHT 4 0 2 1 5 0 11 27 36 11 4
## 60612 60613 60614 60615 60616 60617 60618 60619 60620 60621 60622
## API 1 1 1 0 4 1 0 1 2 0 1
## BLK 245 37 24 154 63 351 15 296 339 339 31
## HISPANIC 23 21 10 2 8 80 60 2 5 2 37
## WHT 24 24 12 4 8 13 22 1 5 2 21
## 60623 60624 60625 60626 60628 60629 60630 60631 60632 60633 60634
## API 1 4 6 2 4 3 8 0 1 0 0
## BLK 426 899 16 110 414 184 17 3 19 1 9
## HISPANIC 134 43 50 44 8 130 29 2 165 2 32
## WHT 21 71 32 27 9 29 64 9 28 0 23
## 60636 60637 60638 60639 60640 60641 60642 60643 60644 60645 60646
## API 3 2 2 0 5 0 0 0 1 3 1
## BLK 389 299 9 37 111 13 13 86 626 38 1
## HISPANIC 12 1 12 107 23 51 10 1 23 17 3
## WHT 6 5 18 17 25 25 5 9 25 13 4
## 60647 60649 60651 60652 60653 60654 60655 60656 60657 60659 60660
## API 1 4 0 0 1 0 0 1 0 8 5
## BLK 19 237 320 32 203 18 2 0 46 24 32
## HISPANIC 87 3 95 11 2 0 4 1 16 31 14
## WHT 8 5 41 10 4 7 7 5 19 23 7
## 60661 60707 60827
## API 1 0 0
## BLK 3 3 31
## HISPANIC 0 7 1
## WHT 4 2 0Preliminary visualizations
Stop Rate by Race
total_pop <- Model_Data %>% filter(Year==2010) %>% select(total_population)
population_by_ethnic <- Model_Data %>%
filter(Year==2010) %>%
mutate_at(c("percent_asian", "percent_black", "percent_hispanic", "percent_white"), function(x){x* total_pop$total_population/100}) %>%
mutate(zip = as.numeric(zip)) %>%
select(percent_asian, percent_black, percent_hispanic, percent_white, zip)
population_by_ethnic## percent_asian percent_black percent_hispanic percent_white zip
## 1 1751.76 1070.52 583.92 6131.16 60601
## 2 321.86 29.26 175.56 907.06 60602
## 3 264.60 35.28 44.10 538.02 60603
## 4 66.40 0.00 12.45 332.00 60604
## 5 4245.24 4494.96 1248.60 14483.76 60605
## 6 362.57 55.78 251.01 2091.75 60606
## 7 4775.46 4021.44 2010.72 13321.02 60607
## 8 7879.40 12607.04 45700.52 11819.10 60608
## 9 2510.92 17576.44 30758.77 10671.41 60609
## 10 2711.80 6585.80 2324.40 26343.20 60610
## 11 4458.30 1188.88 1188.88 22291.50 60611
## 12 1442.08 21991.72 5047.28 7210.40 60612
## 13 2972.16 2972.16 4953.60 37647.36 60613
## 14 3418.15 2734.52 4785.41 56057.66 60614
## 15 2852.15 24447.00 2037.25 10186.25 60615
## 16 19301.49 11877.84 5444.01 11877.84 60616
## 17 0.00 45054.90 30870.95 5840.45 60617
## 18 3896.52 1948.26 45784.11 42861.72 60618
## 19 0.00 62974.34 649.22 649.22 60619
## 20 0.00 71055.41 732.53 732.53 60620
## 21 0.00 32263.68 336.08 336.08 60621
## 22 1628.16 4341.76 15196.16 32020.48 60622
## 23 0.00 28953.72 53649.54 1703.16 60623
## 24 0.00 37318.85 785.66 785.66 60624
## 25 11055.24 3158.64 30007.08 33165.72 60625
## 26 3574.27 12765.25 11744.03 21445.62 60626
## 27 0.00 67407.40 2151.30 1434.20 60628
## 28 1138.33 26181.59 76268.11 10244.97 60629
## 29 6367.46 578.86 15629.22 34152.74 60630
## 30 577.98 288.99 3467.88 24275.16 60631
## 31 4542.25 1816.90 75401.35 9084.50 60632
## 32 134.74 3368.50 5659.08 4311.68 60633
## 33 2908.72 727.18 25451.30 42903.62 60634
## 34 0.00 37522.92 1197.54 399.18 60636
## 35 1930.80 37167.90 965.40 7723.20 60637
## 36 560.33 1680.99 22973.53 30257.82 60638
## 37 0.00 13790.10 70789.18 6435.38 60639
## 38 8367.45 11585.70 8367.45 34113.45 60640
## 39 2829.20 1414.60 38194.20 26877.40 60641
## 40 915.60 2014.32 3845.52 10804.08 60642
## 41 0.00 37321.92 1555.08 11922.28 60643
## 42 0.00 44773.14 1428.93 952.62 60644
## 43 7036.95 7506.08 10320.86 20641.72 60645
## 44 2831.40 283.14 3397.68 21518.64 60646
## 45 2614.38 5228.76 45315.92 33115.48 60647
## 46 0.00 42717.36 908.88 908.88 60649
## 47 0.00 38524.50 20179.50 1834.50 60651
## 48 0.00 20408.16 14880.95 6802.72 60652
## 49 310.86 28288.26 621.72 1243.44 60653
## 50 1588.50 794.25 1429.65 11596.05 60654
## 51 281.73 1972.11 1690.38 24228.78 60655
## 52 2197.44 274.68 3296.16 21150.36 60656
## 53 4766.16 2042.64 4766.16 55151.28 60657
## 54 12690.24 1982.85 7931.40 15862.80 60659
## 55 5467.80 6309.00 7570.80 21450.60 60660
## 56 1635.92 446.16 446.16 4759.04 60661
## 57 1281.36 2989.84 12386.48 24772.96 60707
## 58 0.00 26263.24 1141.88 856.41 60827print(colSums(population_by_ethnic))## percent_asian percent_black percent_hispanic percent_white
## 152464.1 885195.5 785951.5 893236.1
## zip
## 3516807.0data_stops = data.frame(race=c("Asian Pacific Islander", "Black", "Hispanic", "White"), stop_rate=(rowSums(Stops)/colSums(population_by_ethnic %>% select(-zip))))
rownames(data_stops) <- c()
head(data_stops %>% arrange(desc(stop_rate)))## race stop_rate
## 1 Black 0.0080750524
## 2 Hispanic 0.0023182092
## 3 White 0.0009448790
## 4 Asian Pacific Islander 0.0005443903Blacks and Hispanics are stopped at higher rate than Whites and Asians given their population.
Plot of Traffic Stops by Population
l <- population_by_ethnic %>%
mutate(BLK=percent_black, HISPANIC=percent_hispanic, API=percent_asian, WHT=percent_white) %>%
select(BLK, HISPANIC, API, WHT, zip) %>%
gather(key=RACE_CODE_CD, value=population, - zip)
data <- Police.r %>%
filter(ZIP_CD %in% valid_zips) %>%
group_by(ZIP_CD, RACE_CODE_CD) %>%
dplyr::summarise(total.stop = n()) %>%
inner_join(l, by=c("ZIP_CD" = 'zip', 'RACE_CODE_CD')) %>%
mutate(race=RACE_CODE_CD)
ggplot(data, aes(x=population, y=total.stop, colour=RACE_CODE_CD)) + geom_point() + geom_smooth() +
labs(x="Population Count", y="Number of Traffic Stops", title="Traffic Stops By Population", subtitle="Visualization of number of traffic stops as the population of each ethnic group increases for each zipcode") +
guides(colour=guide_legend("Race")) +
theme(plot.title = element_text(hjust=0.5)) ## `geom_smooth()` using method = 'loess'
Plot of Traffic Stops by Percentage of Population
ggplot(data, aes(x=population, y=total.stop, colour=RACE_CODE_CD)) + geom_point() + geom_smooth() +
labs(x="Percent of Population", y="Number of Traffic Stops", fill="Race", title="Traffic Stops By Percentage of Population") +
guides(colour=guide_legend("Race")) +
theme(plot.title = element_text(hjust=0.5))## `geom_smooth()` using method = 'loess'
Police Model 0
Variables Used
Data is only from January 2016 - February 2017 (13/12)
i: ethnic group
j: zipcode locations in Chicago
y[i, j]: observed number of arrests of ethnic group i in area j
n[i, j]: the number of traffic stops in area j for ethnic group i
\(\mu\): grand mean of arrests
\(\alpha\)[i]: the random effect for each ethnic group i
Data
\[y_{ij} \sim Pois(\frac{13}{12} n_{ij} \theta_{ij})\] \[log(\theta_{ij}) = \mu + \alpha_{i} \]
Priors
\[\mu \sim N(0, 0.000001)\] \[\alpha_{i} \sim N(0, 0.000001)\]
log(0.28603)## [1] -1.251659log(0.28603+0.05911) - log(0.28603)## [1] 0.1878534log(0.28603-0.05712) - log(0.28603)## [1] -0.2227678log(0.28603-0.12381) - log(0.28603)## [1] -0.5671433police_model_0 <- "model{
#Data
for (i in 1:4) {
for(j in 1:58){
y[i,j] ~ dpois((13/12) * n[i,j] * theta[i,j])
log(theta[i,j]) <- mu + alpha[i]
}
}
#Priors
mu ~ dnorm(0, 1.0E-06)
for (i in 1:4){
alpha[i] ~ dnorm(0, 0.005)
}
}"
#initial value set
inits = list(mu=-1.251659, alpha=c(0, 0.1878534, -0.2227678, -0.5671433))
police_jags_0 <- jags.model(textConnection(police_model_0), data=list(y=Arrests, n=Stops),
inits=inits)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 232
## Unobserved stochastic nodes: 5
## Total graph size: 597
##
## Initializing model#simulate a sample from the posterior
#note that we specify both mu and tau variables
police_sim_0 <- coda.samples(police_jags_0, variable.names=c("mu", "alpha"), n.iter=40000)
#store the samples in a data frame:
police_data_0 <- data.frame(police_sim_0[[1]])Police Model 1
#estimate initial values
data <- Police.r %>%
filter(ZIP_CD %in% valid_zips) %>%
group_by(ZIP_CD, RACE_CODE_CD) %>%
dplyr::summarise(total.arrest = as.integer(sum(ENFORCEMENT_ACTION_TAKEN_I)), total.stop = n())
head(data)## # A tibble: 6 x 4
## # Groups: ZIP_CD [3]
## ZIP_CD RACE_CODE_CD total.arrest total.stop
## <int> <chr> <int> <int>
## 1 60601 BLK 7 15
## 2 60601 HISPANIC 3 4
## 3 60601 WHT 2 4
## 4 60602 BLK 2 4
## 5 60603 BLK 2 6
## 6 60603 HISPANIC 1 1freq_model <- lm(total.arrest ~ RACE_CODE_CD * total.stop, data = data)
summary(freq_model)##
## Call:
## lm(formula = total.arrest ~ RACE_CODE_CD * total.stop, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.097 -1.148 -0.147 1.273 46.686
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.08564 1.60786 -0.053 0.95758
## RACE_CODE_CDBLK -4.88013 1.84214 -2.649 0.00875 **
## RACE_CODE_CDHISPANIC 0.54742 1.85234 0.296 0.76791
## RACE_CODE_CDWHT 1.12422 1.95807 0.574 0.56655
## total.stop 0.28603 0.49111 0.582 0.56098
## RACE_CODE_CDBLK:total.stop 0.05911 0.49113 0.120 0.90433
## RACE_CODE_CDHISPANIC:total.stop -0.05712 0.49135 -0.116 0.90757
## RACE_CODE_CDWHT:total.stop -0.12381 0.49393 -0.251 0.80235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.55 on 189 degrees of freedom
## Multiple R-squared: 0.9789, Adjusted R-squared: 0.9781
## F-statistic: 1250 on 7 and 189 DF, p-value: < 2.2e-16Variables Used
Data is only from January 2016 - February 2017 (13/12)
i: ethnic group
j: zipcode locations in Chicago
y[i, j]: observed number of arrests of ethnic group i in area j
n[i, j]: the number of traffic stops in area j for ethnic group i
\(\mu\): grand mean of arrests
\(\alpha\)[i]: the random effect for each ethnic group i
\(\beta\)[j]: the random effect for each zipcode location j
\(\epsilon\)[i, j]: the random error for each ethnic group i and location j
Data
\[y_{ij} \sim Pois(\frac{13}{12} n_{ij} \theta_{ij})\] \[log(\theta_{ij}) = \mu + \alpha_{i} + \beta_{j} + \epsilon_{ij}\]
Priors
\[\beta_{0} \sim N(0, 0.000001)\] \[\alpha_{i} \sim N(0, 0.000001)\] \[\beta_{j} \sim N(0, b^2)\] \[\epsilon_{ij} \sim N(0, e^2)\]
Hyperpriors
\[b \sim Gamma(0.5, 0.0005)\] \[e \sim Gamma(0.5, 0.0005)\]
police_model_1 <- "model{
#Data
for (i in 1:4) {
for(j in 1:58){
y[i,j] ~ dpois((13/12)*n[i,j]*theta[i,j])
log(theta[i,j]) <- mu + alpha[i] + beta[j] + epsilon[i,j]
}
}
#Priors
mu ~ dnorm(0,1.0E-06)
for (i in 1:4){
alpha[i] ~ dnorm(0, 0.005)
}
for (j in 1:58){
beta[j] ~ dnorm(0, 1/b^2)
}
for (i in 1:4){
for (j in 1:58){
epsilon[i,j] ~ dnorm(0, 1/e^2)
}
}
b ~ dgamma(0.5, 0.0005)
e ~ dgamma(0.5, 0.0005)
}"
#initial value set
inits = list(mu=-1.251659, alpha=c(0, 0.1878534, -0.2227678, -0.5671433), tau_beta=170, tau_epsilon=170, beta=rep(0.03, 58), epsilon=matrix(rep(0.01, 58*4), ncol=58))
police_jags_1 <- jags.model(textConnection(police_model_1), data=list(y=Arrests, n=Stops),
inits=inits)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 232
## Unobserved stochastic nodes: 297
## Total graph size: 1471
##
## Initializing model#simulate a sample from the posterior
#note that we specify both mu and tau variables
police_sim_1 <- coda.samples(police_jags_1, variable.names=c("mu", "alpha", "beta", "sd_beta", "sd_epsilon", "epsilon"), n.iter=40000)
#store the samples in a data frame:
police_data_1 <- data.frame(police_sim_1[[1]])Running mean plot of beta.1., mu, and alpha.2.
running_mean_plot(police_data_1$beta.1.)
running_mean_plot(police_data_1$mu)
running_mean_plot(police_data_1$alpha.2.)
police_viz <- police_data_1 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "beta")) %>%
filter(grepl("beta", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
mutate(mu = mean(police_data_1$mu)) %>%
mutate(rate = exp(mu + beta)) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
police_model_beta <- left_join(zip_df, police_viz, by="id")
head(police_model_beta)## long lat order hole piece id group zipcode beta
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 0.05094085
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 0.05094085
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 0.05094085
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 0.05094085
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 0.05094085
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 0.05094085
## mu rate zip
## 1 -0.2179287 0.8462099 60601
## 2 -0.2179287 0.8462099 60601
## 3 -0.2179287 0.8462099 60601
## 4 -0.2179287 0.8462099 60601
## 5 -0.2179287 0.8462099 60601
## 6 -0.2179287 0.8462099 60601Map of rate at which traffic stops result in arrests by zipcode (Using Model 3)
model_1_map_beta <- ggplot()
model_1_map_beta <- model_1_map_beta + geom_polygon(data = police_model_beta, aes(x=long, y=lat, group=group, fill=rate), color = "black", size=0.2)
model_1_map_beta <- model_1_map_beta + coord_map()
model_1_map_beta <- model_1_map_beta + labs(x="", y="", title=paste("Rate at Which Traffic Stops Result in Arrests by Zipcode"), fill="total")
model_1_map_beta <- model_1_map_beta + cleanup
print(model_1_map_beta)
Rate of traffic stops by ethnic group
police_viz1 <- police_data_1 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "alpha")) %>%
filter(grepl("alpha", zipcode)) %>%
mutate(mu = mean(police_data_1$mu)) %>%
mutate(rate = exp(mu + alpha))
head(police_viz1)## zipcode alpha mu rate
## 1 alpha.1. -1.255402 -0.2179287 0.2291610
## 2 alpha.2. -1.097255 -0.2179287 0.2684249
## 3 alpha.3. -1.279364 -0.2179287 0.2237350
## 4 alpha.4. -1.344437 -0.2179287 0.2096395